class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 4: Procesamiento masivo de datos geoespaciales ###Procesamiento masivo de datos: Paralelización de rásters José A. Lastra<br> <a href="http://github.com/JoseLastra"> Github: JoseLastra</a><br> <a href="mailto:jose.lastra@pucv.cl"> jose.lastra@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Noviembre 2022</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## Contenidos .pull-left[ - Paralelización: * LibrerÃa parallel * LibrerÃa foreach y doParallel - LibrerÃa terra: * app() * sapp() * lapply() ] .pull-right[ <img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/> © Allison Horst ] --- ## Recordando algunos conceptos -- - Problemas de cálculo ([Jones, 2017](https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html)): * **CPU**: proceso toma mucho tiempo de CPU. * **Memoria**: Demanda demasiada memoria. * **I/O**: toma mucho tiempo la lectura/escritura desde el disco. * **Network**: La red toma demasiado tiempo en transferir. -- - La paralelización viene a "solucionar" los problemas asociados a la CPU, dada la presencia de múltiples núcleos en los pc's modernos. --- ## Paralelización: Conceptos <center><img src="data:image/png;base64,#sock_fork.png" width="900px"/></center> <center>Fuente: Benito, 2021</center> --- ## foreach y doParallel: Ejemplo práctico con rásters -- - Crearemos una lista de archivos **MOD13Q1** correspondientes al Ãndice NDVI desde el año 2000 hasta el 2021. -- - Ejecutaremos dos tareas simples: (1) cortar y (2) enmascarar datos a un área de estudio. -- - Usaremos los siguientes archivos: (1) **san_javier.gpkg**, (2) **MOD13Q1_dates_ts492.csv** y (3) **MOD13Q1_bands (carpeta)** ```r ## Listando archivos bandas <- list.files(path = 'MOD13Q1_bands/',pattern = '*.tif',full.names = T) ## archivo de fechas tabla <- read_csv('MOD13Q1_dates_ts492.csv') ## shape de San Javier shp <- read_sf('san_javier.gpkg') ``` --- ## foreach y doParallel: Ejemplo práctico con rásters -- - Preparación de clúster de procesamiento. ```r ## nombres de salida archivos nombres <- paste0('San_Javier_MOD13Q1/', tabla$dates,'_', tabla$ID,'_San_Javier.tif') ## Cluster set up nCluster <- detectCores(logical = F)-1 #Número de núcleos fÃsicos cl <- makeCluster(nCluster,type = 'PSOCK') #Configuración del cluster registerDoParallel(cl) #registrar para doParallel ``` --- ## foreach y doParallel: Ejemplo práctico con raster -- - Configuración de foreach y detención del clúster ```r foreach (i= 1:length(nombres)) %dopar% { #importar librerÃas a los nodos library(raster) library(magrittr) r <- bandas[i] %>% raster() #crear banda en R m <- r %>% crop(shp) %>% mask(shp) # #cortar y enmascarar #guardar resultado writeRaster(m, filename = nombres[i], overwrite=T, datatype = 'INT2S') } stopCluster(cl) ``` --- <center><img src="data:image/png;base64,#cpu.png" width="950px"/></center> <center>LabGRS, 2022</center> --- ## terra y la paralelización -- - Los ejemplos de paralelización vistos anteriormente, donde se envÃa información a cada nodo, no son compatibles con los objetos de la librerÃa **terra** ([ver más, pág. 6](https://cran.r-project.org/web/packages/terra/terra.pdf)) -- - terra provee funciones que permiten realizar operaciones empleando múltiples núcleos. Siendo **app()** la más común y simple de emplear. ```r ndvi_sj <- list.files(path = 'San_Javier_MOD13Q1/', pattern = '*.tif', full.names = T) %>% rast() ``` --- ## app() -- - Esta función nos permite operar en un **SpatRaster** multi o monobanda para aplicar una función normalmente de resumen (sum, min, max, median, etc.). También podemos operar en **SpatRasterDataset**. -- - **Importante**: dado la eficiencia computacional de **terra** muchas funciones resumen pueden ser aplicadas directamente con un funcionamiento óptimo. -- - Sin embargo, para archivos muy grandes puede ser útil la aplicación de **app()** ```r mean_noApp <- ndvi_sj %>% mean(na.rm = T) #simple mean_app <- app(ndvi_sj, fun = mean, na.rm = T, cores = 5) #5 cores ``` ``` ## Unit: seconds ## expr min lq mean median uq max neval ## simple_ver 1.407792 1.451339 1.450418 1.456992 1.466379 1.469589 5 ## cores_ver 1.411071 1.447485 1.449573 1.450149 1.450531 1.488627 5 ``` -- - **app()** opera los archivos de la misma forma que **apply()**: donde cada banda es una columna en la matriz y los pÃxeles son las filas. --- ## sapp() -- - Si lo que queremos es aplicar una función a cada banda dentro de nuestro **SpatRaster**, podemos emplear *sapp()*. ```r ndvi <- list.files(path = 'MOD13Q1_bands/', pattern = '*.tif', full.names = T) %>% rast() shp_v <- shp %>% vect() #sf to SpatVector #versión simple ndvi_sj <- ndvi %>% crop(shp_v) %>% mask(shp_v) ``` ``` ## Unit: seconds ## expr min lq mean ## ndvi_sj <- ndvi %>% crop(shp_v) %>% mask(shp_v) 1.975383 2.015378 2.057485 ## median uq max neval ## 2.048154 2.100203 2.161207 30 ``` --- ## sapp() ```r # versión con sapp() ## Creando la función crop_mask <- function(x, ...){ r <- x %>% crop(shp_v) %>% mask(shp_v) r } ##aplicación ndvi_sj_sapp <- sapp(ndvi_sj, fun = crop_mask) ``` ``` ## Unit: seconds ## expr min lq mean ## ndvi_sj_sapp <- sapp(ndvi_sj, fun = crop_mask) 6.20839 6.350794 6.419038 ## median uq max neval ## 6.386043 6.478633 6.671329 5 ``` --- ## Observaciones -- - **sapp()** invoca a **sapply()** por detrás y aún no se ha implementado la paralelización para esta alternativa (**Información actualizada al 14 de Octubre de 2022**). -- - Tal cuál lo indica la documentación: *"(...) La paralelización es soportada, pero a menudo no ayuda, e incluso puede ser más lenta."* -- - Justificaciones para paralelizar con terra: * Disponemos de más de 8 núcleos para procesar * O tenemos una función muy compleja (lenta) para aplicar sobre nuestros archivos. -- - Muchas veces el uso de **lapply()** puede ser suficiente. --- ## lapply() con terra -- - Tomemos nuestra lista de archivos originales y la función creada anteriormente para cortar y enmascarar. ```r ndvi_list <- list.files(path = 'MOD13Q1_bands/', pattern = '*.tif', full.names = T) #adaptando la función para lapply crop_mask <- function(x, ...){ x <- rast(x) #SpatRaster r <- x %>% crop(shp_v) %>% mask(shp_v) #aplicación de operación rm(x) # liberar memoria r #entregar objeto } ndvi_crop <- lapply(ndvi_list,FUN = crop_mask) ``` --- ## BibliografÃa complementaria - Bosco, J. (2018). "R Para Principantes". [Ver](https://bookdown.org/jboscomendoza/r-principiantes4/) - Cao, R. y Fernández, R. (2022). "Introducción al procesamiento en paralelo en R". En: "Técnicas de Remuestreo". [Ver](https://rubenfcasal.github.io/book_remuestreo/intro-hpc.html) - Gillespie, C., & Lovelace, R. (2021). "Efficient R programming: a practical guide to smarter programming." O'Reilly Media, Inc. [Ver](https://csgillespie.github.io/efficientR/) - Hijmans, R. J., Bivand, R., Forner, K., Ooms, J., Pebesma, E., & Sumner, M. D. (2022). Package ‘terra’. [Ver](https://brieger.esalq.usp.br/CRAN/web/packages/terra/terra.pdf) - Jones, M. (2017). "Quick Intro to Parallel Computing in R". [Ver](https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html) - McCallum, E., & Weston, S. (2011). "Parallel R". O'Reilly Media, Inc. - Orellana, J. (2018). "HPC con R para investigadores". [Ver](https://bookdown.org/content/1498/) --- class: inverse middle 